load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

location_f="/Mypath/"

fc   = addfile(location_f + "ECEarth_new_20210204_Plumb_250_control_djf.nc", "r")
fw   = addfile(location_f + "ECEarth_new_20210204_Plumb_250_lowice_djf.nc", "r")

;*******************************************;
;Read in variables
;*******************************************
lat=fc->lat
lon=fc->lon
nlat  = dimsizes(lat)
nlon  = dimsizes(lon)

epcx = fc->Fx
epcy = fc->Fy

epwx = fw->Fx
epwy = fw->Fy

;************************
; calculate the anomaly:
;************************

  epx_ano = epwx
  epx_ano = epwx-epcx

  epy_ano = epwy
  epy_ano = epwy-epcy

printVarSummary(epx_ano)

;epy_ano@lat =fc->lat(::8,::16)
;epx_ano@lat =fc->lat(::8,::16)
;-------------------------
; making the plots
;-------------------------
wks = gsn_open_wks("pdf","epflux250hpa_ano_djf_ecearth")

gsn_define_colormap(wks,"BlueDarkOrange18")
plot = new(1,graphic)                          ; create a plot array

  vecres                  = True            ; vector only resources
  vecres@gsnDraw          = False           ; don't draw
  vecres@gsnFrame         = False           ; don't advance frame
  vecres@vcGlyphStyle     = "CurlyVector"   ; curly vectors
  vecres@vcRefLengthF     = 0.02           ; define length of vec ref
  vecres@gsnRightString   = " "             ; turn off right string
  vecres@gsnLeftString    = " "             ; turn off left string
  vecres@tiXAxisString    = " "             ; turn off axis label
  vecres@vcRefAnnoOrthogonalPosF = -.07    ; move ref vector into plot
  vecres@vcLineArrowThicknessF = 2.5
  vecres@vcLabelsOn = False

  vecres@mpCenterLonF     = 290.
  vecres@gsnPolar   = "NH"
  vecres@mpMinLatF            = 0


 vecres@vcRefMagnitudeF  = 1.

;remove very small arrows and the very high latitude ones for nicer plot
 do i=0,nlat-1
 do j=0,nlon-1

 if ( ((abs(epx_ano(i,j)).lt.0.2).and.(abs(epy_ano(i,j)).lt.0.2)) .or. (lat(i).gt.70) ) then
 epx_ano(i,j)=0.
 epy_ano(i,j)=0.
 end if
 
 end do
 end do


plot(0)   = gsn_csm_vector_map_polar(wks,epx_ano(::12,::12),epy_ano(::12,::12),vecres)
;________________________
; creating the panel:
;________________________
 resP            = True

resP@gsnPanelLabelBar    = False;True                ; add common colorbar

resP@gsnPanelBottom   = 0.05                   ; add space at bottom

resP@gsnMaximize    = True                ; maximize plot
gsn_panel(wks,plot,(/1,1/),resP)

